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Abstract 

In the simulation of protoplanetary disc with a power law density 
profile a disc instability is detected. The instability arises only with a 
power law profile and is affected by power index. Thus the impact of 
initial density profile is large within the employed numerical model. 

1 Introduction 

Density profile is the dependence of disc surface density on radius. As it is 
pointed in [I], at present the true density profiles in protoplanetary discs are 
unknown. Nevertheless in most works density profile is thought to be a power 
law with different index. The index should be -1.0 to correspond the observation 
data [2]. In simulation various indexes in the range from -0.5 to -2.5 are being 
used n, 13], 0. 

MMSN model (Minimal Mass Solar Nebula) was proposed in [3]. The den- 
sity of solid particles in this model was obtained by imagining grinding up the 
planets, distributing their mass smoothly with radius and adding up enough gas 
to make Solar composition. The resulting density profile is the following: 



cr(r) = CTi (—^ ) 

Vl a.u./ 



-1.5 



Here ai — 1700 g/cm? , 1 a.u.= 1.5 x 10^^ m. The ratio of gas and solid 
particles mass in MMSN model is the same as in the Solar System (100:1). 
Unfortunately, this model is scarcely applicable to extrasolar planetary systems 
as it follows from observation of T Tauri stars [5] . Therefore the simulation of 
protoplanetary discs with different initial density profiles is conducted [I] , [5] . 

The aim of the present work is to simulate the evolution of the protoplanetary 
disc with different initial density profiles and to compare the results with the 
known simulations of planet system formation. 
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2 Simulation 



Two density profiles were taken as initial profiles for computational experiments: 
solid body profile (as) and power law profile (crp): 




o-s = o-i\/l- 7^ , r < Rd, 



ap = air", a = — 0.5, — 1.5. 

Here Bd is the disc radius and the value of ai is set for the disc mass to be 
equal to the given value. 

The computational experiments were conducted in size-less variables in order 
to decrease round-off errors. The following quantities were chosen as basic 
characteristic parameters for transition to size-less variables: 

— distance from the Sun to the Earth Rq = 1.5 ■ 10^^ m; 

— mass of the Sun M© = 2 • 10^" kg; 

— gravitational constant G = 6.672 • 10~^^ • m^/kg^. 

Corresponding characteristic values of the particle velocity (Vq), time {to), 
potential ($o) and surface density (ao) are written as 



Vo = y^=30 km/s, 

to = ^ = 5 • 10'^ s = 1/6 year, 
'^o-Vq - —5—, 



= -BT- 

In the following text all the parameters are given in size-less units. 

The ratio of central body mass, gas mass and dustsolid particles mass was 
set according to the MMSN model: central body mass Mq = 10.0, gas mass 
was Mq = 1.0 and solid particles mass Mp — 0.01. Both solid particles and gas 
were given the Keplerian velocity vk (in any point of the disc the centrifugal 
force is equal to gravitational one): 

r dr 

Other parameters: initial disc radius Rd = 2.0, radius of computational domain 
Rm = 9.0. 

The dynamics of the solid particle component of protoplanetary disc is de- 
scribed by the Vlasov-Liouvillc kinetic equation. In the following text dustsolid 
particles will be called simply particles. To consider motion of the gas com- 
ponent the equations of gas dynamics are employed. The gravitational field is 
determined by Poisson equation. 
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If wc employ the coUisionless approximation of the mean self- consistent field, 
then Vlasov-Liouville kinetic equation is written in the following form 

Of _^df 

+ vVf + a— = 0, 
at ov 

where /(t, f, v) is the time-dependent one-particle distribution function along 

Ffr . 

coordinates and velocities, a = — V$ H is the acceleration of unit mass 

m 

particle, Ffr is the friction force between gas and particle components of the 
medium. Gravitational potential $ could be divided into two parts: 

^ = + ^2 

where $i presents the potential of protostar. The second part of of potential 
$2 is determined by the additive distribution of the moving particles and gas. 
$2 satisfies the Poisson equation 

A$2 = 47rGSp 

In the case of a flat disc the bulk density of the mobile media Sp = ppart + Pgas 
is equal to zero {ppart is the particle density Pgas 

is the gas density). At the disc 
with the surface density a there is a shear of the normal derivative of potential. 
This shear gives a boundary condition for the normal derivative of potential $2: 

^ = 2.Ga 
oz 

The equations of gas dynamics take the following form: 

| + V(p.1=0 



-Wp + F 



— + {vV)E=-V (pv) + Q+(^F,vj- VW 



where E = ye + ^ j is the density of gas full energy, e = e (p, T) is the internal 

energy of gas, p = p{p, T) is the pressure of gas, W = VT is the heat flow, Q is 
the increase of energy due to chemical reactions and radiation. F is the external 
force which is defined by the following expression 

F = />V$ - kfr [u - v) 

here kfr is the coefficient of friction between gas and particle components of the 
disc, u is the particle velocity, v is the gas velocity. In the case of the flat disc 
the form of equation remains the same with the only exclusion: bulk density p 
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is replaced with surface density cr. In this paper we shall consider only the flat 
disc model. 

Vlasov-Liouville equation is solved by Particles-in-Cells method. To solve 
the equations of gas dynamics Fluids-in-Cells method is employed. Poisson 
equation is solved by a combination of FFT and SOR methods. A detailed 
description of the code could be found in [5] . 

In the computational experiments the cylindrical coordinate system was 
used, grid size is Nu x x Nz = 300 x 256 x 100. The experiments were 
conducted with the MVS-IOOOM multicomputer of the Siberian Supercomputer 
Centre (32 Alpha21264 processors, 833 MHz). 



3 Results 

The most interesting result is that the disc with the massive central body is 
unstable when the initial density profile satisfies the power law. Instability here 
is the loss of axial symmetry in the central are of the disc, as it is shown in 
figure [1] a group of dens gaseous clumps is formed around the central body. 




-0.5 0.5 



Figure 1: Gas density in the central are of the disc. 



It is important because usually the massive central body suppresses all the 
angular instabilities [7]. Furthermore, this instability arises for the disc with 
power law profile and does not arise for the disc with solid body profile. 

In order to show the development of instability the Fourier analysis of gas 
density was conducted along angular direction: 



cr{r,'P,t)^ ^ Sk{r,t)cos (^k(p\ 
k=i ^^^^ ^ 



1 / 27r 

Sk{r,t) ^ ir^ Sk{r,t)cos I —kLp 



SMAxit) ^ uia-KSkir,t), < r < 1 < fc < A^^ 
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Figure 2: Instability development in the disc with power law profile for various 
indexes a. 

Figure [5] shows maximal harmonic amplitude Smax depending on time with 
various power indexes a. 

It should be noticed that zero harmonic (k = 0) is not displayed in figure 
[5] because this harmonic shows instabilities that do not break axial symmetry. 
One can see from figure [2] that the harmonic amplitude increases sufficiently 
with time. 

Now let us consider the behavior of disc with various indexes a. Figure 
shows that with lesser values of a the instability develops faster and harmonics 
have greater amplitude. This fact correspond the results of jX : in their sim- 
ulations discs with lower a formed planets earlier and the planets had greater 
mass. 

It is also stated in [Ij and [3J that for steeper profiles (lower values of a), 
the terrestrial planets are more massive. This result is supported by the figure 
[3] that shows average grain particle density depending on radius. 

4 Conclusion 

It follows from the computational experiment that within the employed model 
of the protoplanetary disc the impact of the initial density profile is large. The 
power law profile leads to the development of the angular instability, while the 
disc with solid body profile remains stable. Moreover, the instability develops 
faster with lower values of power index. The decrease of the power index also 
leads to increase of amplitude of unstable harmonics and to higher mass of gas 
and grain clumps. 
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Figure 3: Gas density in the central are of the disc. 
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